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Abstract 

Quasiparticle bands of the two-dimensional Hubbard model are calculated 
using the Roth two-pole approximation to the one particle Green's function. 
Excellent agreement is obtained with recent Monte Carlo calculations, in- 
cluding an anomalous volume of the Fermi surface near half-filling, which 
can possibly be explained in terms of a breakdown of Fermi liquid theory. 
The calculated bands are very flat around the (tt, 0) points of the Brillouin 
zone in agreement with photoemission measurements of cuprate superconduc- 
tors. With doping there is a shift in spectral weight from the upper band to 
the lower band. The Roth method is extended to deal with superconduc- 
tivity within a four-pole approximation allowing electron-hole mixing. It is 
shown that triplet p-wave pairing never occurs. Singlet d x 2_ y 2-w&ve pairing is 
strongly favoured and optimal doping occurs when the van Hove singularity, 
corresponding to the flat band part, lies at the Fermi level. Nearest neigh- 
bour antiferromagnetic correlations play an important role in flattening the 
bands near the Fermi level and in favouring superconductivity. However the 
mechanism for superconductivity is a local one, in contrast to spin fluctuation 
exchange models. For reasonable values of the hopping parameter the tran- 
sition temperature T c is in the range 10-100K. The optimum doping 5 C lies 
between 0.14 and 0.25, depending on the ratio U/t. The gap equation has a 
BCS-like form and 2A max /kT c ~ 4. 

71.10,+x, 71.27.+a, 74.20.Mn, 74.72.-h 
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I. INTRODUCTION 



In 1957 Bardeen, Cooper and SchriefferB provided, at a stroke, the essentials of a com- 
plete theory of the phenomenon of superconductivity as it was known then and for the next 
thirty years. However since the discovery of high temperature superconductivity in cuprate 
materials! it has been generally, although not universally, believed that a new mechanism 
is operating in these systems. It was established early on that one key element of the BCS 
mechanism, electron pairing, remains.^ However there is mounting evidence that the sym- 
metry of the pairs is d^^taO rather than s-wave and this suggests an electronic mechanism 
for pairing rather than the original BCS phonon-mediated one. 

The common element in all the cuprate superconductors is the Cu0 2 plane in which the 
Cu atoms form a square lattice with an O atom at the midpoint of each pair of nearest- 
neighbour Cu atoms. In the simplest case of La 2 CuO± = (LaO)2Cu02, with the assumption 
of La 3+ and 2 ~ ions, the Cu charge is 2+ corresponding to a 3d 9 configuration. In the 
presence of crystal field splitting one expects doubly occupied d xy , d yz , d zx and d 3z 2_ r 2 
orbitals and a singly-occupied d x i_ y i orbital. In the absence of interaction between the 
electrons the undoped system La 2 CuO± would therefore be a metal whereas it is observed 
to be an antiferromagnetic insulator. This demonstrates the importance of strong repulsive 
Coulomb interaction on the Cu site which tends to localize the d electrons and produces a 
Mott insulator. In the doped system La 2 - x Sr x Cu0 4: , where La 3+ ions are replaced by Sr 2+ , 
the nominal occupation of the x 2 — y 2 Cu orbital is reduced to 1 — x and with x ~ 0.15 the 
system is metallic and superconducting with T c ~ 35K. 

Anderson! was the first to propose that the essence of high temperature superconductiv- 
ity is contained in the two dimensional (2D) square lattice Hubbard mode!0 with repulsive 
on-site interaction U. The atomic orbital in the model may be regarded as a Cu d x 2_ y 2 
orbital hybridized with O p x and p y orbitals on neighbouring sites. 00 In his recent work 
Andersonll! attributes the existence of superconductivity to electron transfer between Cu0 2 
planes. He proposes that electrons in the Cu0 2 plane separate into uncharged spinons and 
charged holons, as they do in the one dimensional Hubbard model, where the electrons 
form a Luttinger liquid. This spin-charge separation inhibits the transfer of unpaired elec- 
trons between planes, but in Anderson's theory this constraint is removed upon pairing. 
The resultant decrease in total kinetic energy favours pairing and drives superconductivity. 
A difficulty with Anderson's theory is the absence of convincing evidence for spin-charge 
separation in two dimensions. 

From the experimental point of view it is true that the critical temperature T c tends 
to increase with the number of adjacent Cu0 2 planes, although Tl 2 Ba 2 CuOg with rather 
isolated Cu0 2 planes has T c = 85K.Ej Some interplanar interaction is presumably essential 
to stabilize superconductivity at finite T even when the primary mechanism, which sets the 
scale of the resultant T c , operates within a plane. An analogous situation is a quasi-2D 
ferromagnet without magnetic anisotropy in which the scale of the Curie temperature is set 
by in-plane exchange interaction, although weak interplanar exchange is required to stabilize 
ferro magnetism at finite T. 

An in-plane mechanism which is advocated by Pines and co-workers,i Scalapino,El and 
MoriyaS is based on electron pairing due to exchange of antiferromagnetic spin fluctuations. 
This type interaction leads inevitably to d x 2_ y 2 pairing and much of the experimental evi- 
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dence cited in its favour is evidence for d x 2_ y 2 pairing rather than for the spin fluctuation 
mechanism itself. In this paper we show how d x 2_ y 2 pairing arises in the 2D Hubbard model 
in a way which is not obviously related to antiferromagnetic spin fluctuations, although 
antiferromagnetic correlations between neighbouring sites are clearly present. The spin fluc- 
tuation mechanism has been strongly criticized by Andersonjil who emphasizes the heuristic 
nature of Pines's model. 

The published ab initio Monte Carlo calculations on the 2D Hubbard model have not 
produced any evidence of superconductivity. El However for the small finite size systems 
considered the electron density cannot be varied continuously and it is possible that the 
favourable case of optimal doping is missed. Very recently Husslein et alM have addressed 
this question by implementing projector Monte Carlo calculations for the 2D Hubbard model 
with nearest and next-nearest neighbour hopping parameters t, t' . This work was motivated 
by the 'van Hove scenario'^ in which, given a suitable interaction between electrons, a high 
T c may emerge from the large density of states at the van Hove singularity associated with 
saddle-points in the band structure. To model a particular system t'/t is tuned so that 
for optimal doping, i.e. with a convenient electron density close to that with the highest 
observed T c , the singularity occurs at the Fermi energy. In earlier wor k§ on the 'van Hove 
scenario' explicitly attractive electron interactions, such as the phonon mediated one, were 
considered and consequently the pairing was s-wave. The remarkable result of the new Monte 
Carlo calculations^ is that d x 2_ y 2 pairing emerges in the repulsive U Hubbard model, but 
only when, for a given t'/t, the doping level is tuned to within ±5% of the optimum one. 
It is clear why superconductivity was missed in earlier Monte Carlo calculations. The new 
Monte Carlo results provide very satisfying confirmation of some of the results reported in 
this paper, which are obtained by a more analytical Green's function method. We will show 
how a combination of on-site electron repulsion and the band structure saddle-point leads 
to d x 2_ y 2 pairing and a fairly high T c . Our approach gives a unified theory of quasiparticles 
in the normal and superconducting states. A slight surprise about the work of Husslein et 
a/.E] is that their modest interaction strength U — 2t leads to sufficient correlation to give 
the effect. 

We use a Green's function decoupling scheme originally due to Linderberg and OhrnH 
and first applied to calculations on the Hubbard model by Roth.^1 The formalism is reviewed 
in Sec. [TT] and in Sec. [TT1] it is applied to the normal paramagnetic state of the 2D Hubbard 
model. It is shown that the Roth two-pole approximation gives excellent agreement with 
the quasiparticle dispersion curves found in recent Quantum Monte Carlo results by Bulut 
et alEj. These authors have shown that their results are consistent with recent angular 
resolved photoemission measurements of the hole-doped cuprates. A natural extension of 
the method to superconductivity, now using a four-pole approximation, is made in Sec. [TV. 
We find a superconducting state in which the gap is determined by a non-standard correlation 
function. In Sec. [V] a gap equation is derived in two ways, yielding a lower and upper bound 
to the gap. It is shown that triplet pairing cannot occur, but that singlet d x 2_ y 2 pairing is 
strongly favoured, with T c = 10 — 100K. The critical temperature T c is strongly dependent 
on doping, with the optimum doping corresponding to the case where the saddle-point of 
the band structure is situated exactly at the Fermi level. This links our work to the 'van 
Hove scenario', although none of its ideas were imposed beforehand, unlike work by other 
authors.B Our results arise from the Hubbard model and the decoupling scheme only. 
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II. THE FORMALISM 



In this paper we wish to determine Green's functions of the Hubbard mo del0 with 
hamiltonian 



h = t E c l* c j* + U J2 n iv n i-v - a* E 



n, 



(2.1) 



<i,j> CT 



including the chemical potential in a standard notation. The equation of motion for a 
retarded Green's function ((A:B)) takes the form 



u((A:B 



A,B 



A,H 



(2.2) 



where < ... > denotes the thermal average. Note the introduction of a new Green's function 
(/ A, H ;-S)/> which implies that an infinite set of equations needs to be solved. A very 

satisfying way to decouple this set of equations is to introduce a set of operators {A n }, which 
are believed to be the most relevant to describe the one particle excitations of the system of 
interest. Formally this assumption is that 



A n} H 



(2.3) 



Here, for the application of the formalism on the 2D Hubbard model in the thermody- 
namic limit for several parameter ratio's U/t, we follow Rothi in choosing two operators of 
Bloch type 



.4 



lkc 



'ha 



1 



e ik.Ri 



(2.4) 



^2ka — ^ka — e^'^ni_ a c ia (2.5) 

to describe the normal, spatially uniform, state of the system. Here L is the number of lattice 
sites in the system. The restriction to spatial uniformity allows discussion of ferromagnetic 
states,ilEil but additional operators are necessary to describe antiferromagnetic statesBcl 
To describe superconducting states, we shall also introduce additional 'hole' operators to 
supplement the 'electron' operators (|2.4j),(|2.5|): 



^ = ^ = ^E^4_ ff , (2.6) 



The coefficients K nm in Eq. ( |2.3j ) are determined by anticommuting both sides of Eq. 
fl2.3| ) with each element of the operator set {A n } and then taking the thermal average. This 
can be written in matrix notation as 



4 



E = KN, 



(2.8) 



where the energy and normalization matrices, E and N, are given by 

E 



A n , H 



4t 

J m 



(2.9) 



N„ 



A 4t 



Combining Eqs. Q, (p|), and flTgP- (|2~l0| ) we obtain 



Eg, 



where G is given by 

G = N(wN-E)" 1 . 

If B = we obtain for the Green's function matrix, whose elements are given by 

G nm (uj) = 

a solution in terms of the matrices E and N: 

GM = N(luN - E) _1 N. 



A ■ A^ 

■"■ni m 



(2.10) 



(2-11) 



(2-12) 



(2-13) 



(2.14) 



This decoupling procedure was first proposed by Linderberg and 6hrn.ll Soon after RothH 
applied this procedure to study ferromagnetism in the 3D infinite U Hubbard model in the 
thermodynamic limit. It can be shown that the formalism is essentially equivalent to the 
Mori-Zwanzig projection technique^ and strongly related to moment methods!! based 
on the assumption that the spectral function is a finite sum of weighted delta functionsBa-B 
The matrix elements E nm and N nm involve correlation functions which should be deter- 
mined self-consistently from the calculated Green's functions as far as possible. Sometimes 
however, as discussed in Sec. [TTI|, one must introduce further Green's functions. The standard 
relationship between correlation function and Green's function may be written as 



BA 



1 

2tt7 



f(u)((A;B)) dco 



(2.15) 



where f(uj) is the Fermi function f(u>) = {er^ + 1) _1 and the contour encircles the real 
axis without enclosing any poles of f(co). The chemical potential \i is determined by the 
condition 



n 



'k a C kcj 



(2-16) 



where n = (n ia ) is the average site occupation per spin. 
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III. THE NORMAL STATE 



For the normal paramagnetic state of the Hubbard model, within the 2-operator set 
{A,r , A 2 t }, the energy and normalization matrices E 2 , N 2 are given byi! 



e^ + Un- fi 



fi + U)n 



(eg - fi + U)n (U - /i)n + epi 2 + n(l - n)Wg 



(3.1) 



N 5 



1 n 

n n 



(3.2) 



The Green's function matrix is readily found using Eq. ( |2.14[ ). Here eg is the unperturbed 
band energy 



<j>i 



and Wr may be written as: 



with 



n(l - n)Wjj = w + , 



W = - J! t (4x C .?>(l - ™i-<7 

<j>i 



n 



]- a > ' 



(3.3) 



(3.4) 



(3.5) 



W\ — \ ((NjNi) — (Nj) (Ni)) 



+ (S j .S, 



(3.6) 



Here Nj = nj a +nj_ a and Sj are the total number operator and the spin operator respectively 
for site j. Sites j and % are nearest neighbours as indicated by the j summation < j >j in 
Eqs. ( ^.3| ) and ( p.5[ ). For the 2D Hubbard model on a square lattice the unperturbed band 
energy is given by 



= 2t (cos k x a + cos fc^a) 



(3.7) 



with a the lattice constant. To model the Cu0 2 plane we take t < so that the bottom 
of the band is at the T-point k x = k y = 0. The occupation number n and the correlation 
functions in wq may be calculated from Green's functions Gu and G\2 by means of Eq. 
( [2.15] ). The correlation functions in w% cannot be determined directly in this way. A natural 
way to determine the density and spin correlation functions in w\ would be to extend the 
formalism to deal consistently with the corresponding two-particle Green's functions. This 
would have the advantage of yielding the spin dynamics of the system and work along 
these lines is in progress. Here, however, we follow Roth's original procedure and introduce 
extra operators Correlations of the form (^AiBi^ can then be calculated within the 
decoupling (|2.3|) by using Eqs. (|2.11|) and ( [2.15|) . The details are given in Roth's paper .0 
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Thus to investigate the normal state of the 2D Hubbard model we follow precisely the 
method Rothi originally applied to the three dimensional cubic lattice. The use of two 
operators, A\ and A 2 corresponds to a two-pole approximation to the Green's function. 
The resulting Green's function therefore consists of two terms, each corresponding to a 
quasiparticle band: 



G llS {u>) = + (3.8) 



with 



1 U(l - 2n) - et + Wt , , 

«ia? = 2 + — 2X~ 5 a 2k = 1 ~ a ^ ( 3 - 9 ) 



o o~ ' ^2fc - ?iig + A fc> l d - iU J 



and 



X S = - e s + W s ) 2 + 4nC% - W n ). (3.11) 



It has been shown that this Green's function conserves the first four moments of the spectral 
density.il The quasiparticle bands uj = £ x £ and u> = £ 2 £ satisfy det (uN 2 — E 2 ) = 0. They 
are plotted in Fig. [1], for U — 8|t| and different site occupations (N) = 2n, along symmetry 
lines in the two dimensional Brillouin zone. The quasiparticle energies are relative to the 
chemical potential \i and comparison is made with the non-interacting band er, relative to 
the non-interacting chemical potential, and with very recent Quantum Monte Carlo results 
by Bulut et a/.@. The agreement with the Monte Carlo results is remarkable and this is 
a strong indication that the Roth two-pole approximation represents the Green's function 
rather well. This is reasonable since the spectral functions of Bulut et a/.il are dominated 
by two peaks which define the quasiparticle energies for each k. Near k = (it, tt) there is also 
an additional broad bump at lower energy, but its weight is very small. One striking feature 
is the flattening at the top of the lower band around k = (ir,ir), extending to the saddle- 
point at k — (vr, 0) and halfway to k — (0,0). Very flat bands have also been observed 
near the (tt, 0) point in recent angular resolved photoemission (ARPES) experiments on 



hole doped cuprate superconductors sue 



. ....... Las Bi 2 Sr 2 CaCu 2 8+s (Bi 2212)S Bi 2 Sr 2 CuO^ 

(Bi2201)Ji YBa 2 Cu 4 8 (YBCO 124), @@ and YBa 2 Cu 3 7 (YBCO 123) SB It is clear 
from our calculations that this feature is associated to a large extent with nearest neighbour 
antiferromagnetic correlations (^Si.Sj^ < 0. Because (^Si.Sj^ is by far the most important 
contributor to Wi for the occupations shown, Wi will be negative, so that W% decreases 
with increasing er. This leads to a narrowing and, near the top, a flattening of the lower 

band. (Si.Sj) increases with occupation, so the effect is pronounced when approaching half- 
filling. The flattening near the top is also responsible for a clear gap between the upper and 
lower band, even for U no larger than the unperturbed bandwidth. The densities of states 
calculated from the one-particle Green's function are shown in Fig. |[ 
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One feature that the Quantum Monte Carlo results appears to share with our results 
is an anomalous Fermi surface volume for (N) > 0.8. In Fig. |3] we plot the Fermi surfaces 
associated with our calculated bands and the non-interacting bands of Fig. [1|. It is clear that 
the volume, or area in two dimensions, of the calculated Fermi surface is always larger than 
the non-interacting one. Since the Quantum Monte Carlo bands agree well with ours along 
symmetry lines, the corresponding Fermi surfaces may be expected to be similar. Certainly 
for (N) = 0.94 it seems that the Quantum Monte Carlo data are definitely incompatible 
with a Fermi surface of normal volume. However, according to Luttinger's theorem0 the 
volume enclosed by the Fermi surface of an interacting Fermi liquid is equal to that of the 
non-interacting system. We would not want to conclude that there is a breakdown of Fermi 
liquid theory on the basis of our own calculations, because, as will be discussed below, the 
two-pole approximation will always lead to a Fermi surface enclosing an enlarged volume 
when only the lower band is occupied. But the Monte Carlo method in principle yields 
exact results for the model, so the apparent anomalous Fermi surface volume must be taken 
seriously. It must be remembered however that the Quantum Monte Carlo calculations^ 
were executed at a finite temperature, T = 0.5 \t\ in this case. A possibility is that the 
enlarged Fermi surface volume is a temperature effect and that the volume shrinks to its 
normal volume as T — > 0, although this means that the form of the Fermi surface and the 
bands would change considerably in the process. In order that the chemical potential moves 
in the right direction with increasing temperature, the Fermi level at T = would have to 
lie above the van Hove singularity where the density of states has a negative gradient. Thus 
electron interaction would have produced a strong distortion of the Fermi surface which, 
although conserving the enclosed volume, changes the topology. Another explanation might 
be that a transition from a Fermi liquid to a non-Fermi liquid, of the type predicted by 
Edwards and Hertz,il is occurring as (N) increases through a critical value N c . In the 
Edwards-Hertz theory, based on a modified Hubbard alloy analogy, a T = phase transition 
occurs at a value of N c < 1 when U is larger than the bandwidth. For (N) < N c we have 
a normal Fermi liquid with a Fermi surface of normal volume defined by a finite Migdal 
discontinuity associated with quasiparticles of infinite lifetime. For (N) > N c there is no 
Migdal discontinuity since all quasiparticles have a finite lifetime, attributed by Edwards 
and Hertz to very strong spin disorder scattering. However clear quasiparticle peaks in the 
spectral functions still exist, just as in the Monte Carlo data, and the points in k space where 
they cross the chemical potential define a 'Fermi surface'. It is foundEl that for (N) > N c 
the volume of this Fermi surface is enlarged compared with the non-interacting one and 
spectral functions and quasiparticle bands have been calculated for the 2D Hubbard model 
within the Edwards-Hertz approach so that a detailed comparison with the Monte Carlo 
data can be made. It is clear on general grounds that, if the 'Fermi surface' volume is 
indeed anomalous, one would not expect a Migdal discontinuity with sharp quasiparticles 
because this is a vestige of the non-interacting situation and would only occur at a surface 
of normal volume. 

We now return to the present Roth approach which cannot describe a sharp transition 
at a critical N c but which gives a smooth crossover from Fermi surfaces of almost normal 
volume at smaller (N) to enlarged surfaces at (N) closer to 1. The major defect of the 
two-hole approximation is that the quasiparticles, whose energies appear to be given rather 
accurately, are not subject to any lifetime broadening. The origin of the enlarged Fermi 
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surface within this approximation is seen most clearly in the limit of very large U where the 
effect is most pronounced. For U — > oo the one-electron Green's function Gn is easily found 
from Eqs. ( |3.8j )-( pTTT| ) to take the form 



G 



I — n 
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uj - (1 - n)eg - nW^ + 
n 

+ u-ne l :-(l-n)W f :-U + fi' (3 ' 12) 

For (N) = 2n < 1 only the lower band, associated with the pole of the first term, is occupied. 
Since each state in this band has weight 1 — n, whereas for a non-interacting case each state 
has weight 1, the Fermi surface volume is increased relative to the non- interacting one by 
a factor (1 — n)^ 1 . For example if (N) = |, the enlargement factor is | and the Fermi 
surface volume corresponds to a half-filled ((N) = 1) non-interacting situation. Since Wt 
depends on k only through eg , as is seen from Eq. fl3.4| ), it follows from Eq. ( ^.10| ) that 



a quasiparticle constant energy surface enclosing a given volume coincides with the non- 
interacting energy surface enclosing the same volume. Thus, with U — > oo, for (N) = | we 
have the square Fermi surface eg = 0, but for finite U, with a less pronounced enlargement 
factor, the square Fermi surface occurs for larger (N). The square Fermi surface plays an 
important role in the superconducting state described in the following sections because, for 
given U, optimum doping (highest T c ) occurs with the occupation (N) corresponding to the 
square Fermi surface. This corresponds to having a van Hove singularity at the Fermi level 
for some finite doping and we have argued above that this must be the case at T = to be 
consistent with the high temperature Monte Carlo data. Thus even if the anomalous Fermi 
surface volume of the Roth method is spurious we believe the topology of the Fermi surface 
is described quite well. 

We wish to emphasize again the remarkable agreement shown in Fig.l between the 
quasiparticle bands calculated by the Roth two-pole approximation, using precisely her 
method of evaluating various correlation functions which occur, and those given by the 
Monte Carlo calculations of Bulut et Our calculations agree with the Monte Carlo 

data not only for the location of the bands, but also for the spectral weight distribution over 
the two bands. This is shown in Fig. |] where we show, for several points along the k x = k y 
line in the Brillouin zone, our two delta functions with weights and a 2 g. These weights 
compare well with the areas under the peaks of the spectral functions found by Bulut et 
al£j While the spectral weight distribution is dispersionless for U — > oo, for intermediate U 
(U = 8\t\) we find a strong dispersion. We also find that there is a strong weight transfer 
upon doping from the upper band u = to the lower band u = ^g, as shown in Fig. |^. This 
weight transfer is particularly strong around the (ir, tt) point for intermediate U. A strong 
weight transfer from high energy scale to low energy scale upon doping has been observed 
in electron loss experiments!-^ and optical spectroscopy experimentsea on high temperature 
cuprate superconductors. Applying strong-coupling perturbation theory, Eskes et a/.0 find 
a similar weight transfer. 

We may mention here that the Roth procedure, again followed precisely, also gives a good 
description of the occurence of ferromagnetism in the 2D Hubbard model. The ferromagnetic 
region of the ((N),t/U) plane is found to bell somewhat larger than the best variational 
estimates,^ but the result is again remarkably good considering the subtlety of the energy 



9 



balance in this problem and the simplicity of the method. This contrasts strongly with the 
Hubbard I approach,^ which corresponds to setting = 0. RothH pointed out that the 
constant shift Wq in W% [see Eq. ( |3.4| )1 for the minority spin band is essential to stabilize 
the state of complete spin allignment for large U near half-filling. Recently Hewson and 
WassermanS showed how an approach equivalent to the Roth one leads to the results of 
slave-boson mean field theory in the Anderson impurity model; here again the shift in the 
impurity arising from wq is essential to place the resonant state correctly near the Fermi 
level. These successes of the present approach encourage us to consider another subtle 
problem, superconductivity, in the next section. We shall see that an anomalous correlation 
function which characterizes the superconducting state is generated automatically, just like 
the correlation functions in Wt. 

IV. THE SUPERCONDUCTING STATE 

To discuss superconductivity we need to mix electron and hole operators and evaluate 
anomalous correlation functions in which particle number is not conserved. In this context 
we note that the procedure described in Sec. [II], applied on the BCS reduced Hamiltonian 
and using the operator set {c^ ff , } yields the whole BCS formalismlll in Nambu-Gor'kov 
form.E! As indicated in Sec. 0, for the Hubbard model we use the four operators ( |2.4| )- ([2.7|) 
and hence obtain a four-pole approximation to the Green's functions. The 4x4 energy and 
normalization matrices E4 and N4 may each be partitioned into four 2x2 matrices. The 
upper left 2x2 matrices are just E 2 and N 2 , given by Eqs. (|3~l"l ) and (|3.2|) , and the lower 
right 2x2 matrices are easily shown to be — E 2 and +N 2 . The elements of the off-diagonal 
blocks involve anomalous correlation functions and follows: 

A 13 = A 24 = A 31 = A 42 = 0, (4.1) 

N u = ~N 23 = N* 41 = -N* 2 = (c^ew) , (4.2) 

E 13 = E* 31 = U (a-vdr) , (4.3) 

E u = E; 2 = (eg - n) (c laCi - a ) , (4.4) 

E23 = #41 = ( U - 6 k ~ I 1 ) ( C i-° C i°) 

+ t ( C i-crCla + Cl„ a C ia ) , (4.5) 

<l>i 



<l>i 

+ t (nicrCi-orCla ~ 71^0^0^) , (4.6) 

<l>i 
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E i2 = E t 



Jk.{Ri-R t ) 



t J 



n io C l-o C lo + n l-crCi-a C 



<l>i 



(4.7) 



These expressions simplify considerably when we introduce the symmetry associated with 
either singlet or triplet pairing. The anomalous correlation functions are matrix elements 
between states which differ by the addition of one pair and the spin function associated with 
a pair is symmetric or antisymmetric for triplet or singlet pairing respectively. It follows 
that if all spin labels in the correlation function are reversed, i.e. a — > —a, the sign changes 
in the singlet case but not in the triplet case. 

We first consider triplet pairing. Owing to the spin symmetry the second terms in E24 
and E42 vanish and E13 = En = E23 = 0. Similarly all the elements of the off-diagonal 
blocks in N4 vanish. Furthermore, denoting the correlation function in the first term of E24 
by fin, we have 

fill \H>i—<jClaCl—a ^lo^ia^i—a) 
(^crQ— ctQcj ~t~ ^l — a^i — cr^icr) 

= -Ai- (4-8) 



Using this symmetry we write: 



fin 



fi for Ri-Ri = (a, 0) 
-fi for Ri-Ri = (-a, 0) 
±0 for Ri-R l = (0, a) 
Tfi for Ri-Ri = (0,-a). 



(4.9) 



Thus 



where 



Similarly 



E- 



21 



ft? = 2t/5(sin k x a ± sin k y a) 



E 



12 



(4.10) 



(4.11) 



(4.12) 



It will be seen shortly that /3g is essentially the gap function and the nodes at k x = ±k y 
indicate p-wave symmetry. This is as expected for triplet pairing. It is very interesting that 
the gap is determined not by a simple pair wave function (ci- a ci a ), but by the more subtle 
correlation function fin. This quantity cannot be evaluated directly from our Green's func- 
tions. The simplest way to evaluate it, and hence obtain a gap equation, is by a factorization 
procedure described in Sec. |V| for the case of singlet ci-wave pairing. We shall not pursue 
the p-wave case further here, because this procedure, and others we have investigated, does 
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not yield a non-zero solution for the gap. The reason for this will be pointed out in Sec. [V]. 
We therefore find that p-wave superconductivity does not occur in the 2D Hubbard model. 

We now turn to singlet pairing and discuss particularly the case of <i-wave symmetry. 
We defer a discussion of s-wave symmetry to another occasion, partly because we believe 
it is less likely and partly because the possibility of on-site pairing (ci- a Ci a ) means that all 
elements E nm in the off-diagonal blocks are non-zero. The <i-wave case is simpler because the 
(i-symmetry of the pair wave function ensures that (cj-o-Qo-) = and that the summations of 
the pair wave function (ci^ a ci a ) and of {n ia Ci_ a ci a ) over sites I, which are nearest neighbours 
of site i, give zero. The consistency of these assertions can be checked by expressing these 
correlation functions in terms of the resultant Green's function. Thus E ri = E u = E 2 s = 
and, as in the p-wave case, all elements of the off-diagonal blocks in N4 vanish. In a similar 
way the second terms in £"24 and £42 vanish and we are left with just the first terms of £"24 
and £42, as in the p-wave case. Denoting the correlation function in the first term of E 2 4 
by Ju, we find, introducing the sign changes on spin reversal in the singlet case, that the 
equation analogous to Eq. ( |4.8| ) is: 

lil = (ni-aClaCl-a + TL^C^Ci-a) = 7 W . (4.13) 

Using the ci-wave symmetry we write 



la 



7 for Ri — Ri— (±a, 0) 
-7 for Ri - Ri = (0, ±a) 



(4.14) 



Hence 



E 2 4 = 7g, #42 = 7, 



(4.15) 



where 



ik.iRt-Ri) 



with 



In = g(cos k x a — cos k y a) 



(4.16) 



9 = 2t 7 - 



(4.17) 



We call 7£ the gap-function and g the gap-function amplitude. The matrices E 4 and N 4 
now take the partitioned form 



E2 



il 





7S 



(4.18) 



N 2 
N, 



(4.19) 



where E 2 and N 2 are given by Eqs. (|3.1|) and ( |3.2| ). We shall have a superconducting 
state if 7g 7^ and the central problem we shall presently address is to find a gap equation 
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satisfied by 7^ and search for a non-zero solution. For the moment we assume 7^ 7^ and 
describe how to calculate the Green's function matrix and discuss the general nature of the 
quasiparticle bands. 

The Green's function matrix must be calculated from Eq. (|2.14 ). Here we concentrate 
on the element G13 = ((c£ CT ; c_£ )) from which we can determine the pair wavefunction 
(cj-o-Qo-). The poles of the Green's function give the quasiparticle bands. It is convenient to 
partition the matrix (cjN 4 — E4) into four 2x2 matrices: 



Then 



wN 4 - E 4 



A B 
C D 



(4.20) 



G(u) = N 4 (cuN 4 - E 4 ) _1 N 4 



N 2 
N 2 



(A-BD- l C) 1 —A~ X B (D — CA~ l B) 1 
-D- x C{A-BD- x G)~ l (D-CA^B)' 1 



N 2 
N 2 



(4.21) 



and the calculation of G13 only involves the upper right-hand block —A l B{D — CA 1 B) . 
It is now straightforward to show that 



G 13 H = -1%(E 12 - nE 11 ) 2 /D{u), 



D(u) = det (N 4 u - E 4 ) 



(4.22) 



where 



'u - E n )(um - E 22 ) - (oon - E 12 ) 
x ({u + E u )(um + E 22 ) - (oon + E 12 )' 

En, E\2 and £"22 may be substituted from Eq. ( |3.1| ) and we find 

G n (u) = - ll :(n(l-n)U) 2 /D(u J ). 



(4.23) 



(4.24) 



The quasiparticle bands in the superconducting state satisfy D(u) = and we denote them 
by oj = ±E 1 j: J u = ±E 2 %, with E 2 ^ > > 0. Thus 



D{u) = n\l - nf{u 2 - E 2 ,){u 2 - E\A 



(4.25) 



When |7g|/n(l - n) < - we have 



E 



lk 



e 2 - + 



12 ' 



2k ^lk 



(4.26) 



E, 



2k 



^2k 



n 2 (l-n) 2 £ 2 r -£ 



12 ' 
'life 



(4.27) 
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and, from Eq. Q3.1J) , En = eg + Un — /i. Here £ 2 g with ^ < £ 2 g are the quasiparticle 
bands in the normal state given by Eq. ( 3.10 ). The form of E 2 ~ is familiar from BCS-theory 

and the superconducting gap at points k on the normal state Fermi surface ^ = is given 
by |Ag| where, combining Eq. ( OB ) and Eq. fl4.26| ), 

g (cos k x a — cos k y a) (eg + Un — //) 



n(l — n) 



(4.28) 



*2k 



— * 

The fc-dependence of Ag is entirely due to the (cos k x a — cos fc y a)-factor, since the second 
factor is constant over the Fermi surface. So, over the Fermi surface, 



Ag = G (cos k x a — cosk y a) 



with 



G 



(4.29) 



(4.30) 



We call G the gap amplitude. Recall that we call g, introduced in Eq. ( [4.171) , the gap- 
function amplitude. G tends to g/(l — n) as U — > oo. Examples of quasiparticle bands 
in the superconducting state are given later in Fig. [| with the gap-function amplitude g 
determined by a method described in Sec. [V]. 

In calculations which relate to the superconducting state it is reasonable to assume that 
Wg, which appears in the matrix E 2 [whose elements appear in Eq. ( (4.24| )], takes the values 
it has in the normal state (where 7g = 0) at T = 0. Thus the effect of superconductivity on 
the quasiparticle bands, and the temperature dependence of these bands, is entirely due to 
the gap function 7g. The correlation function (c_j:_ cr c^ is related to G\z(uS) by the use of 
Eq. so that, using Eq. (ggg ) and Eq. ( ^25[) , 



— k — <T KG 



1 

27Ti 



-IkU 2 



-du. 



Hence 



where 



F(a,b) 



c-s- A) = ~JkU 2 F (E& Ex) , 

tanh (h/3aj tanh (\ftb) 



2 (b 2 - a 2 



At T = (ft — > 00) this becomes 



C -k-a C ka 



-UP 2 



2E lk E 2k ( E lk + E 2k) 

The pair wave function in real space may then be calculated using the relation 



ik.lRi-Ri 



C -k-a C kc 



(4.31) 



(4.32) 



(4.33) 



(4.34) 



(4.35) 
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V. THE GAP EQUATION 



In the (i-wave case the gap is determined by the gap-function 7^ and we need a self- 
consistent way of finding the gap-function amplitude g which appears in Eq. (|4.16 ). As 



pointed out in the previous section the gap is not determined by a simple pair wave function 
(ci-aCio), but by the correlation function ju given by Eq. (|4.13|) , which may be written as 

lil = ((4-a c l-<r + c L c io) Ci-aCla) • (5.1) 

In the d-wave case this may also be written as 

7« = (c\ a Ci a Ci- a C io ^ + (civCiaCi-crClJ) , (5.2) 

which shows explicitly that 7^ = 7^. This is not expressible directly in terms of our Green's 
functions and various possibilities suggest themselves. A satisfactory procedure might be to 
extend the formalism to two-particle Green's functions as suggested in connection with the 
correlation functions appearing in Wi, Eq. ( p.6|) . This remains to be investigated. Alter- 
natively we can follow the procedure we adopted, following Roth, for the latter correlation 
functions and introduce additional Green's functions containing one- and three-particle op- 
erators. This is the method we shall use here, since it is consistent with our treatment of 
the normal state which, by comparison with Monte Carlo results, has been shown to work 
well. However there is one complication. Whereas for the correlation functions appearing in 
the normal paramagnetic state the procedure is unambiguous,!! this is not so for 7^. The 
four-operator correlation function may be related to a Green's function containing a one- 
and three-particle operator in four ways and in the case of ju these are distinct and yield 
different results. It is sufficient to consider one of the correlation functions in Eq. Q5.2|) . 
Thus (cfoCitrCi-vCio^ may be calculated using any of the following Green's functions: 

(a) {(c ilT \ c\ a Ci^ a Ci a yj , (b)((c 4 _ a ; c^c^c^)), 

(c) ((c l( ,] C^Ci^Ci^)) , (d)((i; Q CT Q_ a Q CT )). 

Use of (a) or (b) splits up the product Ci a Ci- a and therefore does not preserve the important 
result that 7^ — > when site double occupation is suppressed as U — > 00. Consequently 
both of these procedures, which yield similar but not identical results, fail at large U and 
are generally expected to overestimate the gap. This result is closely related to the factor- 
ization (^yCieCi-aCioS = (cj a c ia ) (ci^ a ci a ) and in fact use of (a) or (b) is exactly equivalent 
to this factorization as U ^00. Rather than use (a) or (b) precisely we develop the related 
factorization method in Sec. [V A[ This has the advantage of simplicity and shows clearly 



the structure of the gap equation and why <i-wave pairing occurs and p-wave pairing does 
not. Although it fails for large U, we believe it is useful in providing an upper estimate of 
the gap and T c at finite, intermediate U. 

Use of the Green's functions (c) or (d) preserves the property ju — > as U ^00 and is 
therefore the prefered procedure for large U. At smaller U it seems likely that the correct 
results may lie between those obtained from (c) or (d) and those from (a), (b), or the 
factorization procedure. We therefore regard the gap and T c obtained using (c) or (d) as a 
lower estimate of those quantities. It is remarkable that although analyses using (c) and (d) 
are formally different the numerical results for the gap are almost identical. In Sec. |VB| we 
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will therefore give details of the gap equation for case (c) only. Of the four possibilities (c) 
has the advantage of giving the correct U — > oo limit for the correlation function and of a 
more balanced distribution of creation and destruction operators than (d). 

There is a reassuring consistency in that all the procedures discussed yield <i-wave pairing, 
not p-wave pairing, and in all cases the gap and T c as function of site occupation (N) are 
sharply peaked at an optimum doping corresponding to the square Fermi surface = 0. 
These results are presented in Sees. |V A| and |V B| for the two basic methods which provide 
upper and lower estimates of the magnitude of the gap and T c . 



A. The factorization procedure 



The factorization procedure approximates ju, given by Eq. (|5.1|), as 



(5.3) 



The symmetry j u = j ti is retained but, as discussed above, the products c io -Cj_ CT and c^c^ 
are split up so that in general 7^ does not tend to zero when site double occupation is 
suppressed as U — > 00. However the averages (ci- a ci a ) are always zero due to <i-wave 
symmetry, so that double occupation within a pair is suppressed, but not between pairs. 
In the special case of half-filling ((N) — 1), ju — > as U ^00, because the correlations 
cl_ a ci-. a ) and (c] a c ia \ tend to zero. We may rewrite Eq. ( |5.3| ) as 



7iz = 2m (ci- a ci a ) 



where 



ni = {ci_ a ci- a 
It is straightforward to show that 



C la C i<J 



Hit 



k 



(5.4) 



(5.5) 



(5.6) 



where = ( c \ a c ka) an d z is the number of nearest neighbours, i.e. z = 4 for the 2D square 
lattice. 

We are now in a position to find the equation which determines the gap-function ampli- 
tude g in Eq. ( |4.16| ) self-consistently. From Eqs. ( |4 . 1 6| ) and ( |5.4| ) we get 



(5.7) 



By using Eqs. ( |4.32| ) and ( |4.35| ) we find: 

2 ni tU 2 



7fe 



J] cos [k. (Mi - ^)] E Tg-cos [q. (r\ - R t )} F (E l? , E 2tf ) 



(5.8) 



(0, 
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Clearly the quasiparticle energies Ei$, E 2 q are symmetric in q x and q y , because they satisfy 
the equation D(u)=0, where D(uS) is defined by Eq. ( 4.23Q and \^q] 2 and the elements of 
E 2 have this symmetry. Now substituting Eq. (|4.16|) , 

7 9 -=5f (cos q x a - cos q y a) , (5.9) 

into Eq. (|5.8|) and exploiting the symmetry between q x and q y , we find: 

2nitU 2 2 
7^ = -1% — 2 — ( cos Q* a ~ cos % a ) E K E ift E W ■ (5.10) 
q 



Thus for a non-zero solution 7^ must satisfy the gap equation 

! 2r2 1 t[/ 2 
L 

k 



(cos k x a - cos k y df F {E^, E 2 ^j , (5.11) 



where 7^ is contained in the quasiparticle energies E^, E 2 %. At T = 0, F(a, b) takes the 
form [2ab(a + b)]" 1 , so that the gap equation becomes: 

mtU 2 ^ (cos k x a — cos k v a) 2 

1 - J2 — t ^r- ( 5 - 12 ) 



L k E lk E 2k { E lk + E 2k 

Since Yj% e k = an< ^ n k ls a m onotonically decreasing function of er, it follows from Eq. 
( |5.6| ) that nit < 0. Also E 2 ^ > E x t > 0, so that a solution with g 7^ may occur. At 
finite temperature T Eq. Q5.ll ) determines g as a function of temperature, and hence the 



critical temperature T c where g = 0. Eq. Q5.12 ) has a rather BCS-like form, particularly 



if we allow ourselves to consider large U where E 2 ^ ~ U and the factor E 2 ^ (E^ + E 2 ^j 
in the denominator in Eq. (|5.12|) cancels with the factor U 2 in the numerator. In this 
approximation the superconductivity is driven by the kinetic energy term —nit, but for large 
U, where the approximation of Sec. [VB| is more appropriate, this is effectively replaced by 
a quantity of order t 2 /U. 



In the case of p-wave pairing, which was briefly considered in Sec. [TV], the equation 
corresponding to Eq. ( p. 12 ) is of the form 



riitU 2 (sin k r a ± sin k„a) 2 
1 = ^— - V- —, V — L x . 5.13 

U E lk E 2k( E lk + E 2k) 

The right-hand side is now negative and no solution is possible. The change of sign from 
the ci- wave case is due to the factor i in Eq. (|4.10D , which does not appear in Eq. ( [4.151 ). 
Thus p-wave superconductivity does not occur. 

As mentioned earlier, in calculating superconducting properties such as the gap ampli- 
tude and its temperature dependence we neglect the unimportant temperature dependence 
of rii and and evaluate these with g = and T = 0. The ci-wave gap amplitude G at 
T = over the Fermi surface calculated from Eqs. ( |5.12| ) and ( |4.30| ) is shown in Fig. [7] 
as a function of occupation (N) for different values of U/\t\. Clearly, from Eq. (|4.28|) the 



maximum gap occurs at points on the Fermi-surface where k x = or k y = 0. For each U/\t\ 
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the gap maximum, corresponding to optimum doping, occurs for the occupation which has 
the square Fermi surface = 0. The flat saddle-point feature of the normal state quasipar- 
ticle bands (Fig. [I]) around (n, 0) then lies precisely at the Fermi level. If the occupation 
is increased (underdoping with holes) the saddle-point lies below the Fermi level and the 
Fermi surface is a closed hole-like one around (tv,tt). For overdoping the saddle-point lies 
above the Fermi level and the Fermi surface is a closed electron-like one around (0, 0). This 
situation is unchanged in the approximation of Sec. |V Jj|, which is more appropriate for large 



U. It has close similarities with the van Hove scenario of Newns et alM and has implications 
for transport properties which we explore in the discussion section. 

Quasiparticle bands in the superconducting state are obtained by solving the equation 
D(u>) = 0, where D(u) is given by Eq. flOp an d the bands are given to a good approxi- 
mation by Eqs. (|4.26|) and ( 4.27|) . Examples of these bands are shown in Fig. []. The upper 



Hubbard bands to = ±E 2 ^ are too far from the Fermi level to appear in this figure. 

In Fig. ^ we print the values of the pair wave function {c%- a ci a ), calculated from Eq. ( 4.85 ) 



for several Ri — Ri on the lattice. It turns out that for optimum doping the pair wave function 
is strongly peaked if i and I are neighbouring sites. Away from optimum doping this peak is 
less pronounced and the pair wave function extends over several lattice sites. This behaviour 
of the pair wave function reflects the local nature of the pairing mechanism. 

The temperature dependence of the gap-function amplitude g(T) is calculated from Eq. 
( |5.11p and in Fig. | results for G(T) are shown for U/\t\ = 12 and various occupations 



including the optimum doping case of (N) = 0.76. It is clear from Fig. ^that 2A max //c B T c ~ 
4 for optimum doping and under-doping, but larger for overdoping. The units of g and k^T 
in Figs. [7] and [| are 2\t\, so that if \t\ takes a reasonable value of 0.5 eV, the units are eV. 
Then T c at optimum doping for U/\t\ = 12 is 125 K. As discussed earlier we regard these 
results for the gap and T c as upper bound estimates. They are strongly reduced in the 
approximation of the next section, which we believe yields lower bound estimates, but the 
qualitative behaviour is unchanged. 

B. A gap equation valid for large U 

As discussed above, we obtain a gap equation by writing the correlation functions which 
appear in Eq. ( |5.2|) as 



[cj a c ia Ci^ a ci a ) = (B K ci a ) , (5.14) 
with 

B\% — CfoCifj-Ci—fj-, (5.15) 
and thus expressing it in terms of the Green's function ((ci a ; Bii)}. From Eq. ( |5.2j ) we have 

7iz = (Biic icT ) + (-B K Q<r) • (5.16) 
It is convenient to introduce Bloch operators by writing 

[Bucu,)=L- 1 ^e^(B tt c s \. (5.17) 
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Also, using Eqs. ( 2.11 ) and ( |2.15 ) and noting that = Cr we have 



where 



I m {k) = — f f(uj)G lm (k,uj)dw. 



It is straightforward to show that 



A - 


Bu 


+> 


A - 


Bu 


+ ) 




Bu 


+ ) 


A - 


Bu 


+ ) 



o. 

L -i/2 e iZA ( ni _ aCi(rCi _ a ) 
L-VV^ ( ni - 



, r-1/2 ifefl, / t J „. 



where n\ is defined by Eq. ( |5.5| ) and, similarly, 



rrii 



C-i&Tli—fjCifj 



(5.18) 
(5.19) 

(5.20) 
(5.21) 
(5.22) 

(5.23) 
(5.24) 



So rrii can be calculated straightforwardly from Gu and it is sufficient to use its value in 
the normal state at T = 0. This may be written as 



(5.25) 



in the case when only the lower quasiparticle band is occupied. It is important to note that 
for large U, niit ~ t 2 /U, whereas riit ~ t. The expressions for Gi m are obtained in a similar 
way to that used to find Gi 3 in Sec. |V[ We find 



Gu 

G13 

Gu 



Un\l - nf {u + fa) (w + ^)/D{u), 
-^Un 2 {l -n)(u + e n -u,+ U)/D(u), 
7fe^n(l - n) (w + eg - /i + Un)/D{u). 



(5.26) 
(5.27) 
(5.28) 



We note also that for <i-wave symmetry the correlation function on the right of Eq. (|5.21 ) 
is equal to (BuCiJ^. It is also important note that, owing to the factors 7^ in Eqs. ( |5.27|) 

and ( |5.28|) , h(k) and I±{k) change sign under interchange of k x and k y . Consequently on 
substituting for the correlation functions on the right of Eq. (|5.18|) from Eqs. (|5.20|) - (|5.23|) 
and substituting the result in Eq. ( |5.17| ), we find that the contribution from the second term 
on the right of Eq. ( |5.23|) vanishes on summing over k. The result is 
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h{k)ni + J 4 (A;)(ni - 



(5.29) 



where 



C 



It follows from Eqs. ( fO§ ) and ( |5TT6| ) that 
2tC 



<9< 



^ COS fc.(-R; - -Ri) X! 7g COS <f.(-R; - i?, 



7g 



(5.30) 



(5.31) 



This is of the same form as Eq. ( |5.8|) and, on using Eq. (|5.9| ) as before, we find the gap 
equation 



2tC 



(cos k x a — cos k y a) 



7fc 



(5.32) 



It is sufficient to evaluate the constant C in the normal state at T = and we find 

-1 



C 



(5.33) 



where 9{x) — 1 if x > and = otherwise. As U — > 00, C — > 1 — n. At T 
( |5.26| )-( f5T2~8| ) and Eq. G5.19| ), we obtain the gap equation in the form 



0, using Eqs. 



1 = 



UtC 



n(l — n)L 



(cos k x a — cos k y ( 

k 



[(1 - n)ni - mi] - e^) + Um\U 

E lk E 2k( E lk + E 2k) 



(5.34) 



This is similar in form to the gap equation, Eq. ( |5.12| ), obtained using the factorization 
method. Both terms in the numerator of the last factor of Eq. ( |5.34|) are of order t, since 
mi ~ t/U. If we were to retain only the Uniin term and take C — 1 — n, its value at 
U = 00, we would obtain the form of Eq. (|5.12|) exactly, but with m\ replacing n\. Thus, 
from the discussion following Eq. ( |5.12| ), superconductivity at large U is driven by a term of 
order t 2 /U. This seen clearly in Fig. 0, where the gap amplitude G over the Fermi surface 
is plotted as a function of occupation (N) for different values of U/\t\. The peak value 
of G, corresponding to optimum doping with the square Fermi surface as before, increases 
rapidly with decreasing U/\t\. At U = 4|t| the peak value is about one tenth of the rather 
constant peak value (see Fig. |?|) found by the factorization procedure. Qualitatively nothing 
is changed and the quasiparticle bands are just like those of Fig. |[ but with a smaller gap. 
The temperature dependence of the gap amplitude g(T) is calculated from Eq. ( |5.32| ) and 
in Fig. [TT] the result for G(T) is shown for U = A\t\ at optimum doping. If \t\ = 0.5eV the 
corresponding T c is 10K. 

As discussed earlier we expect the correct values for the gap and T c to lie somewhere 
between the values obtained from the two approximations used. Thus for the physically 
reasonable parameters \t\ = 0.5eV, U = 2eV, we have shown that ci-wave superconductivity 
occurs with T c in the range 10 — 100K at optimum doping 5 C = 0.14 . 
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VI. CONCLUSIONS 



We have systematically applied Roth's decoupling method^ to the two dimensional Hub- 
bard model, investigating both the normal and superconducting states. Although there are 
intrinsic defects in the method we believe we are able to draw some significant conclusions. 
In the normal state the method corresponds to a two-pole approximation to the one-particle 
Green's function, so that for a given wave-vector k and spin a the spectral function consists 
of two delta-function peaks, the sum of whose weights is unity, corresponding to states in 
the upper and lower Hubbard bands. For site occupation (N) < 1 only the lower band is 
occupied and the absence of any incoherent spectral weight below the Fermi level necessarily 
implies a Fermi surface of volume greater than that appropriate to a non-interacting system 
or Fermi liquid. Nevertheless we have shown in Sec. |TTT| that there is very remarkable agree- 
ment between calculated bands and those obtained in Monte Carlo calculations by Bulut 
et aZ.;0 in particular for U = 8\t\ and several values of (N) the position of the Fermi level, 
and hence the nature the Fermi surface, is in agreement. This implies that for (N) > 0.8 
the anomalous Fermi surface volume is not merely a result of the Roth approximation, but 
is a real effect. This could be caused by the finite temperature in the Monte Carlo calcu- 
lations, but can possibly also be explained on the basis of the Edwards-Hertz theoryffl 
in which for U larger than the bandwidth, a transition from a Fermi liquid to a non-Fermi 
liquid with enlarged Fermi surface volume occurs as (N) increases through a critical value 
N c . However the two-pole approximation omits the quasiparticle lifetime broadening which, 
in the Edwards-Hertz theory occurs even at the Fermi level for (N) > N c . Our calculated 
bands also feature almost dispersionless bands around the (n, 0) points as is observed in re- 
cent ARPES experiments on hole doped cuprate superconductors.11^0 Another feature our 
calculations share with experiments@@ is the weight transfer from the upper to the lower 
band upon doping. For intermediate U this transfer is fc-dependent and has a maximum at 
the (n, 7r) point. 

In Sec. [V] it is shown how c/-wave superconductivity, but not p-wave, occurs. An im- 
portant conclusion is that reasonable transition temperatures, in the range 10 — 100K, only 
occur for U <• 6|i| and in a fairly small range of occupation (N) near optimum doping where 
the saddle-point feature of the quasiparticle band structure at (n, 0) is placed at the Fermi 
level. For U ^ 6|t| the system is expected to be a Fermi liquid for all (N), so that the 
anomalous volume of the Fermi surface is spurious in this regime. However this defect of the 
method is rather fortunate since it enables us to obtain the saddle-point feature of the band 
at the Fermi level for a doped case ((N) < 1) without recourse to next nearest neighbour 
hopping. Also the absence of lifetime broadening in the two-pole approximation is a correct 
feature near the Fermi level in the Fermi liquid regime. We therefore believe our results 
are significant and closely parallel the Monte Carlo results of Husslein et a/S on <i-wave 
superconductivity in the tt' Hubbard model. The association of optimum doping with a 
saddle-point at the Fermi level (see also Ref. |21|) is consistent with the generically different 
behaviour of doped superconducting cuprates in many normal state properties at exactly 
optimum doping, as opposed to those at under- or overdoping. It is well knownlll that a 
saddle-point at the Fermi level leads to marginal Fermi liquid transport properties such as 
a normal state resistivity varying linearly with temperature. A recent discussion of ther- 
mopower0 strengthens the case for the van Hove scenario in the cuprate superconductors. 
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Even though the method of Sec. [VB| still gives superconductivity near optimum doping for 
U > 8\t\, but with very low T c , it is unlikely that the superconducting state would survive 
the introduction of finite quasiparticle lifetimes at the Fermi level which might arise from 
non-Fermi liquid behaviour in this region. 

The present work suggests that high temperature superconductivity arises in moderately 
correlated systems which may be doped so that the Fermi level approaches a van Hove 
singularity in the quasiparticle band structure. Correlation helps by flattening the bands, 
accentuating the singularity, but must not be allowed to reduce site double occupation too 
markedly. Nearest neighbour antiferromagnetic correlations are particularly important both 
in flattening the bands and in promoting double occupation. The correlations important for 
superconductivity in the present model are local ones and the mechanism is different from 
spin-fluctuation-exchange models. 

In future work we propose to explore the possibility of s-wave pairing, particularly its 
degree of anisotropy if it is found to occur. Work is in progress on extending the Roth method 
consistently to two-particle Green's functions, particularly with the aim of calculating the 
dynamical susceptibility to) and investigating the spin dynamics of the system. 
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FIG. 1. The quasiparticle bands along the triangle (0,0) - (vr,7r) - (vr,0) for several occupations 
< N >. The solid lines are the calculated ones for U = 8\t\; the circles are the results from the 
Quantum Monte Carlo calculations by Bulut et oi.0 for the same parameters. The dashed line 
indicates the band in the noninteracting (U = 0) case for the given occupation, (a) < N >= 0.75; 
(b) < N >= 0.87; (c) < N >= 0.94; and (d) < N >= 1.0 (half-filling). 
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FIG. 2. The density of states for U = 8\t\ and the occupations as in figure l(a)-(d). 
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FIG. 4. The single particle spectral weight A(k,E), along the k x = k y line, from our calcula- 
tions. The height of the bars represents the weight of the delta functions. The weights of the delta 
functions compare well with the area's under the peaks of Fig. 1(c) in Bulut et.al, Phys. Rev. B 
50, 7215 (1994). (< TV >= 0.87; U = 8\t\) 




FIG. 5. The spectral weight of the lower band for several occupations < N >. Note that 
the sum of the spectral weights of the two bands is 1, so that this figure shows a weight transfer 
from the upper to the lower band upon doping. (U = 8\t\) 
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FIG. 6. (a) The 2 bands close to the Fermi surface in the superconducting state, calculated 
using the factorization procedure. Note the gap in the neighbourhood of the (ir, 0) point and the 
zero gap on the k x = k y diagonal, reflecting the d-wave symmetry. Here U = 8\t\ and < N >= 0.795 
(nearly optimum doping), (b) shows the neighbourhood of the (ir, 0)-point in detail. The normal 
state electron and hole bands are shown in thin lines. 
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FIG. 7. The gap amplitude G as a function of the occupation < N > for several ratio's U/\t\, 
calculated using the factorization procedure. (\t\ = 0.5eV). 
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FIG. 8. Values of the pair wave function (ci- a ci a ) for several lattice vectors Ri — R[. The 
pair wave function value is printed above the lattice points and the lattice vector is given, as a 
coordinate pair in units of the lattice constant, below the lattice points. U = 12\t\ . (a) < N >= 0.72 
(overdoping) , (b) < N >= 0.76 (optimum doping), and (c) < N >= 0.82 (underdoping) . 
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FIG. 9. The gap amplitude G as a function of temperature for U = 12\t\ = 6eV for several 
occupations < N >, calculated using the factorization procedure. Note that the maximum gap 
A = AG at optimum doping 5 C (=0.24 here). 
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FIG. 10. The gap amplitude G as a function of the occupation < N > for several ratio's U/\t\, 
calculated using a decoupling of the gap-function that is correct in the limit U — > oo. (\t\ = 0.5eV). 
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FIG. 11. The gap amplitude G as a function of temperature for U = 4|t| = 2eV for optimum 
doping 5 C = 0.14, calculated using a decoupling of the gap function that is correct in the limit 
U -> oo. 
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